% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function PV_ldgh(grid,base,lambda,u,q,k,l,outfile,test_name)
% PV_ldgh(grid,base,lambda,u,q,k,l,outfile,test_name)
%
% Write a (lambda,u,q) solution of the LDG-H method in a format
% readable by paraview (.vtk).
%
% To avoid matching the nodes of the hybrid mesh with those of
% the primal and dual variables, the output is divided in two files.
% 1) the file outfile-lambda.vtk contains the hybrid variable,
% represented as a fully discontinuous field;
% 2) the file outfile-uomega.vtk contains the primal and the dual
% variables, structured as discontinuous fields and
% interpolated on the same grid.
%
% See also: interpolate.

 % First the hybrid variable

 % Interpolate lambda to the required degree
 [Xi,Ti,L,cell_type] = interpolate(grid,base,k, ...
        reshape(lambda,[base.nk,grid.ns]),"lambda");
 % add th third component for 2D results
 if(size(Xi,1)==2)
   Xi(3,:,:) = 0;
 end

 % Initialize the output file
 fid = fopen([outfile"-lambda.vtk"],'w');
 
 % Preamble
 line = '# vtk DataFile Version 2.0';
 fprintf(fid,'%s\n',line);
 line = [test_name"-lambda"];
 fprintf(fid,'%s\n',line);
 line = 'ASCII';
 fprintf(fid,'%s\n',line);
 line = 'DATASET UNSTRUCTURED_GRID';
 fprintf(fid,'%s\n\n',line);

 % points
 npoints = size(Xi,2)*size(Xi,3);
 line = ['POINTS ',num2str(npoints),' float'];
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%f ',Xi(:,i,j));
     fprintf(fid,'\n');
   end
 end

 % cells
 line = ['CELLS ',num2str(size(Ti,2)),' ', ...
                  num2str((size(Ti,1)+1)*prod(size(Ti,2)))];
 fprintf(fid,'\n%s\n',line);
 for i=1:size(Ti,2)
   fprintf(fid,'%i',size(Ti,1));
   fprintf(fid,' %i',Ti(:,i)-1);
   fprintf(fid,'\n');
 end

 % cell type
 line = ['CELL_TYPES ',num2str(size(Ti,2))];
 fprintf(fid,'\n%s\n',line);
 fprintf(fid,'%i ',cell_type*ones(1,size(Ti,2)-1));
 fprintf(fid,'%i\n',cell_type); % write the last one followed by \n

 % Write the data
 line = ['POINT_DATA ',num2str(npoints)];
 fprintf(fid,'\n%s\n',line);

 line = 'SCALARS lambda float 1';
 fprintf(fid,'%s\n',line);
 line = 'LOOKUP_TABLE default';
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%f ',L(i,j));
     fprintf(fid,'\n');
   end
 end
 fprintf(fid,'\n');

 fclose(fid);

 % Now the two remaining variables

 % Initialize the output file
 fid = fopen([outfile"-uomega.vtk"],'w');
 
 % Preamble
 line = '# vtk DataFile Version 2.0';
 fprintf(fid,'%s\n',line);
 line = [test_name"-uomega"];
 fprintf(fid,'%s\n',line);
 line = 'ASCII';
 fprintf(fid,'%s\n',line);
 line = 'DATASET UNSTRUCTURED_GRID';
 fprintf(fid,'%s\n\n',line);

 % Interpolate to the required degree
 [Xi,Ti,U,cell_type] = interpolate(grid,base,k, ...
        reshape(u,[base.pk,grid.ne]),"phi");
 [Xi,Ti,O,cell_type] = interpolate(grid,base,k, ...
        reshape(q,[base.mk,grid.ne]),"omega");
 if(size(Xi,1)==2)
   Xi(3,:,:) = 0;
 end
 if(size(O,1)==2)
   O(3,:,:) = 0;
 end

 % points
 npoints = size(Xi,2)*size(Xi,3);
 line = ['POINTS ',num2str(npoints),' float'];
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%f ',Xi(:,i,j));
     fprintf(fid,'\n');
   end
 end

 % cells
 line = ['CELLS ',num2str(size(Ti,2)),' ', ...
                  num2str((size(Ti,1)+1)*prod(size(Ti,2)))];
 fprintf(fid,'\n%s\n',line);
 for i=1:size(Ti,2)
   fprintf(fid,'%i',size(Ti,1));
   fprintf(fid,' %i',Ti(:,i)-1);
   fprintf(fid,'\n');
 end

 % cell type
 line = ['CELL_TYPES ',num2str(size(Ti,2))];
 fprintf(fid,'\n%s\n',line);
 fprintf(fid,'%i ',cell_type*ones(1,size(Ti,2)-1));
 fprintf(fid,'%i\n',cell_type); % write the last one followed by \n

 % Write the data
 line = ['POINT_DATA ',num2str(npoints)];
 fprintf(fid,'\n%s\n',line);

 line = 'SCALARS u float 1';
 fprintf(fid,'%s\n',line);
 line = 'LOOKUP_TABLE default';
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%f ',U(i,j));
     fprintf(fid,'\n');
   end
 end
 fprintf(fid,'\n');

 line = 'VECTORS omega float';
 fprintf(fid,'%s\n',line);
 for j=1:size(Xi,3)
   for i=1:size(Xi,2)
     fprintf(fid,'%f ',O(:,i,j));
     fprintf(fid,'\n');
   end
 end
 fprintf(fid,'\n');

 fclose(fid);
  
return


